function global_dif_Heat = mycode_Heat(flag_rcp, flag_model, flag_period)

%%% Calculate lake heat release feedback changes. The difference between the present-day period and the pre-industrial or future projections

% flag_rcp: The difference between the flag_rcp (pre-industrial or future projections) and present-day period 
% #2 = picontrol        #3 = rcp26      #4 = rcp60      #5 = rcp85

% flag_model: The results of specific lake model or ensemble mean
% #1 = albm     #2 = lake      #3 = simstrat-uog		#4 = vic-lake		#5 = ensemble mean

% flag_period: The seasonal or annual results
% #1 = DJF       #2 = MAM       #3 = JJA      #4 = SON      #5 = annual


swup_rcpMean_ensemble = evalin('base','swup_rcpMean_ensemble');
lwup_rcpMean_ensemble = evalin('base','lwup_rcpMean_ensemble');
SH_rcpMean_ensemble = evalin('base','SH_rcpMean_ensemble');
rcp = evalin('base','rcp');
period = evalin('base','period');

flag_var = 1;   % only 1

abs_rate_swup = 0.23;
abs_rate_lwup = 0.9;
abs_rate_SH = 1;

if flag_model < 5
    
    present_Heat = swup_rcpMean_ensemble.( rcp{ 1 } ).( period{flag_period} ){ flag_model }{ flag_var } * abs_rate_swup + ...
                    lwup_rcpMean_ensemble.( rcp{ 1 } ).( period{flag_period} ){ flag_model }{ flag_var } * abs_rate_lwup + ...
                    SH_rcpMean_ensemble.( rcp{ 1 } ).( period{flag_period} ){ flag_model }{ flag_var };
    future_Heat = swup_rcpMean_ensemble.( rcp{flag_rcp} ).( period{flag_period} ){ flag_model }{ flag_var } * abs_rate_swup + ...
                    lwup_rcpMean_ensemble.( rcp{flag_rcp} ).( period{flag_period} ){ flag_model }{ flag_var } * abs_rate_lwup + ...
                    SH_rcpMean_ensemble.( rcp{flag_rcp} ).( period{flag_period} ){ flag_model }{ flag_var };

else
    
    present_Heat =  cat(3, swup_rcpMean_ensemble.( rcp{ 1 } ).( period{flag_period} ){ 1 }{ flag_var }, ...
                        swup_rcpMean_ensemble.( rcp{ 1 } ).( period{flag_period} ){ 2 }{ flag_var }, ...
                        swup_rcpMean_ensemble.( rcp{ 1 } ).( period{flag_period} ){ 3 }{ flag_var } ) * abs_rate_swup + ...
                    cat(3, lwup_rcpMean_ensemble.( rcp{ 1 } ).( period{flag_period} ){ 1 }{ flag_var }, ...
                        lwup_rcpMean_ensemble.( rcp{ 1 } ).( period{flag_period} ){ 2 }{ flag_var }, ...
                        lwup_rcpMean_ensemble.( rcp{ 1 } ).( period{flag_period} ){ 3 }{ flag_var } ) * abs_rate_lwup + ...
                    cat(3, SH_rcpMean_ensemble.( rcp{ 1 } ).( period{flag_period} ){ 1 }{ flag_var }, ...
                        SH_rcpMean_ensemble.( rcp{ 1 } ).( period{flag_period} ){ 2 }{ flag_var }, ...
                        SH_rcpMean_ensemble.( rcp{ 1 } ).( period{flag_period} ){ 3 }{ flag_var } ) * abs_rate_SH;
    future_Heat =   cat(3, swup_rcpMean_ensemble.( rcp{flag_rcp} ).( period{flag_period} ){ 1 }{ flag_var }, ...
                        swup_rcpMean_ensemble.( rcp{flag_rcp} ).( period{flag_period} ){ 2 }{ flag_var }, ...
                        swup_rcpMean_ensemble.( rcp{flag_rcp} ).( period{flag_period} ){ 3 }{ flag_var } ) * abs_rate_swup + ...
                    cat(3, lwup_rcpMean_ensemble.( rcp{flag_rcp} ).( period{flag_period} ){ 1 }{ flag_var }, ...
                        lwup_rcpMean_ensemble.( rcp{flag_rcp} ).( period{flag_period} ){ 2 }{ flag_var }, ...
                        lwup_rcpMean_ensemble.( rcp{flag_rcp} ).( period{flag_period} ){ 3 }{ flag_var } ) * abs_rate_lwup + ...
                    cat(3, SH_rcpMean_ensemble.( rcp{flag_rcp} ).( period{flag_period} ){ 1 }{ flag_var }, ...
                        SH_rcpMean_ensemble.( rcp{flag_rcp} ).( period{flag_period} ){ 2 }{ flag_var }, ...
                        SH_rcpMean_ensemble.( rcp{flag_rcp} ).( period{flag_period} ){ 3 }{ flag_var } ) * abs_rate_SH;
end

if flag_rcp == 2
    a = future_Heat;
    future_Heat = present_Heat;
    present_Heat = a;
else
end
clear a

global_dif_Heat = future_Heat - present_Heat;      


end
